hud_group_zip <- function(inc_lim, indata){
  
  phia <- c()
  phib <- c()
  theta <- c()
  k <- c()
  fp <- c()
  geoid <- c()
  cbsa <- c()
  tot <- c()
  res <- c()
  
  for(i in 1:nrow(indata)){
    
    a=NA
    b=NA
    phia[i] <- NA
    phia[i] <- NA
    
    geoid[i] <- indata$GEOID[i]
    cbsa[i] <- indata$CBSA[i]
    tot[i] <- indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+ 
      indata$phi5[i]+indata$phi6[i]+indata$phi7[i]+indata$phi8[i]+
      indata$phi9[i]+indata$phi10[i]+indata$phi11[i]+indata$phi12[i]+
      indata$phi13[i]+indata$phi14[i]+indata$phi15[i]+
      indata$phi16[i]
    
    if(inc_lim[i] < 10000){
      fp[i] = indata$phi1[i]
    } else if(inc_lim[i] >= 10000 & 
              inc_lim[i] < 14999.999 & 
              indata$phi1[i] > 0){
      a = 10000
      b = 14999
      phia[i] = indata$phi1[i]
      phib[i] = indata$phi1[i]+ indata$phi2[i]
      if(indata$phi2[i] >= 0.001 & phib[i]!=1){
        theta[i] = (log(1-phia[i])-log(1-phib[i]))/(log(b)-log(a))
        k[i] = ((phib[i] - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi2[i] >= 0.001 & phib[i]==1){
        theta[i] = (log(1-phia[i])-log(1-0.999))/(log(b)-log(a))
        k[i] = ((0.999 - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi2[i] < 0.001){
        fp[i] = phia[i]
      }
    } else if(inc_lim[i] >= 15000 & 
              inc_lim[i] < 19999.999 & 
              indata$phi1[i]+ indata$phi2[i] > 0){
      a = 15000
      b = 19999
      phia[i] = indata$phi1[i]+ indata$phi2[i]
      phib[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]
      if(indata$phi3[i] >= 0.001 & abs(1-phib[i]) > 0.001){
        theta[i] = (log(1-phia[i])-log(1-phib[i]))/(log(b)-log(a))
        k[i] = ((phib[i] - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi3[i] >=0.001 & abs(1-phib[i]) < 0.001){
        theta[i] = (log(1-phia[i])-log(1-0.999))/(log(b)-log(a))
        k[i] = ((phib[i] - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      }else if(indata$phi3[i] < 0.001) {
        fp[i] = phia[i]
      }
    } else if(inc_lim[i] >= 20000 & 
              inc_lim[i] < 24999.999 & 
              indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i] > 0){
      a = 20000
      b = 24999
      phia[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]
      phib[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]
      if(indata$phi4[i] >= 0.001 & abs(1-phib[i]) > 0.001){
        theta[i] = (log(1-phia[i])-log(1-phib[i]))/(log(b)-log(a))
        k[i] = ((phib[i] - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi4[i] >= 0.001 & abs(1-phib[i]) < 0.001){
        theta[i] = (log(1-phia[i])-log(1-0.999))/(log(b)-log(a))
        k[i] = ((0.999 - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi4[i] < 0.001){
        fp[i] = phia[i]
      }
    } else if(inc_lim[i] >= 25000 & 
              inc_lim[i] < 29999.999 & 
              indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ 
              indata$phi4[i] > 0){
      a = 25000
      b = 29999
      phia[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]
      phib[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+ 
        indata$phi5[i]
      if(indata$phi5[i] >=0.001 & abs(1-phib[i]) > 0.001){
        theta[i] = (log(1-phia[i])-log(1-phib[i]))/(log(b)-log(a))
        k[i] = ((phib[i] - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi5[i]>=0.001 & abs(1-phib[i]) < 0.001){
        theta[i] = (log(1-phia[i])-log(1-0.999))/(log(b)-log(a))
        k[i] = ((0.999 - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi5[i] < 0.001){
        fp[i] = phia[i]
      }
    } else if(inc_lim[i] >= 30000 & 
              inc_lim[i] < 34999.999 & indata$phi1[i]+ indata$phi2[i]+ 
              indata$phi3[i]+ indata$phi4[i]+ indata$phi5[i]>0){
      a = 30000
      b = 34999
      phia[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+
        indata$phi5[i]
      phib[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+ 
        indata$phi5[i]+indata$phi6[i]
      if(indata$phi6[i] >=0.001 & abs(1-phib[i]) > 0.001){
        theta[i] = (log(1-phia[i])-log(1-phib[i]))/(log(b)-log(a))
        k[i] = ((phib[i] - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi6[i] >=0.001 & abs(1-phib[i]) < 0.001){
        theta[i] = (log(1-phia[i])-log(1-0.999))/(log(b)-log(a))
        k[i] = ((0.999 - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi6[i]< 0.001){
        fp[i] = phia[i]
      }
    } else if(inc_lim[i] >= 35000 & 
              inc_lim[i] < 39999.999 & indata$phi1[i]+ indata$phi2[i]+ 
              indata$phi3[i]+ indata$phi4[i]+indata$phi5[i]+indata$phi6[i]>0){
      a = 35000
      b = 39999
      phia[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+
        indata$phi5[i]+indata$phi6[i]
      phib[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+ 
        indata$phi5[i]+indata$phi6[i]+indata$phi7[i]
      if(indata$phi7[i] >=0.001 & abs(1-phib[i]) > 0.001){
        theta[i] = (log(1-phia[i])-log(1-phib[i]))/(log(b)-log(a))
        k[i] = ((phib[i] - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi7[i] >=0.001 & abs(1-phib[i]) < 0.001){
        theta[i] = (log(1-phia[i])-log(1-0.999))/(log(b)-log(a))
        k[i] = ((0.999 - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi7[i] < 0.001){
        fp[i] = phia[i]
      }
    } else if(inc_lim[i] >= 40000 & 
              inc_lim[i] < 44999.999 & indata$phi1[i]+ indata$phi2[i]+ 
              indata$phi3[i]+ indata$phi4[i]+indata$phi5[i]+indata$phi6[i]+
              indata$phi7[i]>0){
      a = 40000
      b = 44999
      phia[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+
        indata$phi5[i]+indata$phi6[i]+indata$phi7[i]
      phib[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+ 
        indata$phi5[i]+indata$phi6[i]+indata$phi7[i]+indata$phi8[i]
      if(indata$phi8[i] >= 0.001 & abs(1-phib[i]) > 0.001){
        theta[i] = (log(1-phia[i])-log(1-phib[i]))/(log(b)-log(a))
        k[i] = ((phib[i] - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi8[i] >= 0.001 & abs(1-phib[i]) < 0.001){
        theta[i] = (log(1-phia[i])-log(1-0.999))/(log(b)-log(a))
        k[i] = ((0.999 - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi8[i] < 0.001){
        fp[i] = phia[i]
      }
    } else if(inc_lim[i] >= 45000 & 
              inc_lim[i] < 49999.999 & indata$phi1[i]+ indata$phi2[i]+ 
              indata$phi3[i]+ indata$phi4[i]+indata$phi5[i]+indata$phi6[i]+
              indata$phi7[i]+indata$phi8[i]>0){
      a = 45000
      b = 49999
      phia[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+
        indata$phi5[i]+indata$phi6[i]+indata$phi7[i]+indata$phi8[i]
      phib[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+ 
        indata$phi5[i]+indata$phi6[i]+indata$phi7[i]+indata$phi8[i]+
        indata$phi9[i]
      if(indata$phi9[i]>=0.001 & abs(1-phib[i]) > 0.001){
        theta[i] = (log(1-phia[i])-log(1-phib[i]))/(log(b)-log(a))
        k[i] = ((phib[i] - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi9[i]>=0.001 & abs(1-phib[i]) < 0.001){
        theta[i] = (log(1-phia[i])-log(1-0.999))/(log(b)-log(a))
        k[i] = ((0.999 - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi9[i]< 0.001){
        fp[i] = phia[i]
      }
    } else if(inc_lim[i] >= 50000 & 
              inc_lim[i] < 59999.999 & indata$phi1[i] + indata$phi2[i] + indata$phi3[i] + 
              indata$phi4[i] + indata$phi5[i] + indata$phi6[i] + indata$phi7[i] + 
              indata$phi8[i] + indata$phi9[i]>0){
      a = 50000
      b = 59999
      phia[i] = indata$phi1[i] + indata$phi2[i] + indata$phi3[i] + indata$phi4[i] + 
        indata$phi5[i] + indata$phi6[i] + indata$phi7[i] + indata$phi8[i] + 
        indata$phi9[i]
      phib[i] = indata$phi1[i] + indata$phi2[i] + indata$phi3[i] + indata$phi4[i] + 
        indata$phi5[i] + indata$phi6[i] + indata$phi7[i] + indata$phi8[i] + 
        indata$phi9[i] + indata$phi10[i]
      if(indata$phi10[i]>=0.001 & abs(1-phib[i]) > 0.001){
        theta[i] = (log(1-phia[i])-log(1-phib[i]))/(log(b)-log(a))
        k[i] = ((phib[i] - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi10[i]>=0.001 & abs(1-phib[i]) < 0.001){
        theta[i] = (log(1-phia[i])-log(1-0.999))/(log(b)-log(a))
        k[i] = ((0.999 - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi10[i] < 0.001){
        fp[i] = phia[i]
      }
    } else if(inc_lim[i] >= 60000 & 
              inc_lim[i] < 74999.999 & indata$phi1[i]+ indata$phi2[i]+ 
              indata$phi3[i]+ indata$phi4[i]+indata$phi5[i]+indata$phi6[i]+
              indata$phi7[i]+indata$phi8[i]+indata$phi9[i]+indata$phi10[i]>0){
      a = 60000
      b = 74999
      phia[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+
        indata$phi5[i]+indata$phi6[i]+indata$phi7[i]+indata$phi8[i]+
        indata$phi9[i]+indata$phi10[i]
      phib[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+ 
        indata$phi5[i]+indata$phi6[i]+indata$phi7[i]+indata$phi8[i]+
        indata$phi9[i]+indata$phi10[i]+indata$phi11[i]
      if(indata$phi11[i]>=0.001 & abs(1-phib[i]) > 0.001){
        theta[i] = (log(1-phia[i])-log(1-phib[i]))/(log(b)-log(a))
        k[i] = ((phib[i] - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi11[i]>=0.001 & abs(1-phib[i]) < 0.001){
        theta[i] = (log(1-phia[i])-log(1-0.999))/(log(b)-log(a))
        k[i] = ((0.999 - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi11[i] < 0.001){
        fp[i] = phia[i]
      }
    } else if(inc_lim[i] >= 75000 & 
              inc_lim[i] < 99999.999 & indata$phi1[i]+ indata$phi2[i]+ 
              indata$phi3[i]+ indata$phi4[i]+indata$phi5[i]+indata$phi6[i]+
              indata$phi7[i]+indata$phi8[i]+indata$phi9[i]+indata$phi10[i]+
              indata$phi11[i]>0){
      a = 75000
      b = 99999
      phia[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+
        indata$phi5[i]+indata$phi6[i]+indata$phi7[i]+indata$phi8[i]+
        indata$phi9[i]+indata$phi10[i]+indata$phi11[i]
      phib[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+ 
        indata$phi5[i]+indata$phi6[i]+indata$phi7[i]+indata$phi8[i]+
        indata$phi9[i]+indata$phi10[i]+indata$phi11[i]+indata$phi12[i]
      if(indata$phi12[i]>= 0.001 & abs(1-phib[i]) > 0.001){
        theta[i] = (log(1-phia[i])-log(1-phib[i]))/(log(b)-log(a))
        k[i] = ((phib[i] - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi12[i]>= 0.001 & abs(1-phib[i]) < 0.001){
        theta[i] = (log(1-phia[i])-log(1-0.999))/(log(b)-log(a))
        k[i] = ((0.999 - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi12[i] < 0.001){
        fp[i] = phia[i]
      }
    } else if(inc_lim[i] >= 100000 & 
              inc_lim[i] < 124999.999 & indata$phi1[i]+ indata$phi2[i]+ 
              indata$phi3[i]+ indata$phi4[i]+ indata$phi5[i]+indata$phi6[i]+
              indata$phi7[i]+indata$phi8[i]+indata$phi9[i]+indata$phi10[i]+
              indata$phi11[i]+indata$phi12[i]>0){
      a = 100000
      b = 124999
      phia[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+
        indata$phi5[i]+indata$phi6[i]+indata$phi7[i]+indata$phi8[i]+
        indata$phi9[i]+indata$phi10[i]+indata$phi11[i]+indata$phi12[i]
      phib[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+ 
        indata$phi5[i]+indata$phi6[i]+indata$phi7[i]+indata$phi8[i]+
        indata$phi9[i]+indata$phi10[i]+indata$phi11[i]+indata$phi12[i]+
        indata$phi13[i]
      if(indata$phi13[i] >= 0.001 & abs(1-phib[i]) > 0.001){
        theta[i] = (log(1-phia[i])-log(1-phib[i]))/(log(b)-log(a))
        k[i] = ((phib[i] - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi13[i] >= 0.001 & abs(1-phib[i]) < 0.001){
        theta[i] = (log(1-phia[i])-log(1-0.999))/(log(b)-log(a))
        k[i] = ((0.999 - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi13[i] < 0.001){
        fp[i] = phia[i]
      }
    } else if(inc_lim[i] >= 125000 & 
              inc_lim[i] < 149999.999  & indata$phi1[i]+ indata$phi2[i]+ 
              indata$phi3[i]+ indata$phi4[i]+ indata$phi5[i]+indata$phi6[i]+
              indata$phi7[i]+indata$phi8[i]+indata$phi9[i]+indata$phi10[i]+
              indata$phi11[i]+indata$phi12[i]+indata$phi13[i]>0){
      a = 125000
      b = 149999
      phia[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+
        indata$phi5[i]+indata$phi6[i]+indata$phi7[i]+indata$phi8[i]+
        indata$phi9[i]+indata$phi10[i]+indata$phi11[i]+indata$phi12[i]+
        indata$phi13[i]
      phib[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+ 
        indata$phi5[i]+indata$phi6[i]+indata$phi7[i]+indata$phi8[i]+
        indata$phi9[i]+indata$phi10[i]+indata$phi11[i]+indata$phi12[i]+
        indata$phi13[i]+indata$phi14[i]
      if(indata$phi14[i] >=0.001 & abs(1-phib[i]) > 0.001){
        theta[i] = (log(1-phia[i])-log(1-phib[i]))/(log(b)-log(a))
        k[i] = ((phib[i] - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi14[i] >=0.001 & abs(1-phib[i]) < 0.001){
        theta[i] = (log(1-phia[i])-log(1-0.999))/(log(b)-log(a))
        k[i] = ((0.999 - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi14[i] < 0.001){
        fp[i] = phia[i]
      }
    } else if(inc_lim[i] >= 150000 & 
              inc_lim[i] < 199999.999 & indata$phi1[i]+ indata$phi2[i]+ 
              indata$phi3[i]+ indata$phi4[i]+ indata$phi5[i]+indata$phi6[i]+
              indata$phi7[i]+indata$phi8[i]+indata$phi9[i]+indata$phi10[i]+
              indata$phi11[i]+indata$phi12[i]+indata$phi13[i]+
              indata$phi14[i]>0){
      a = 150000
      b = 199999
      phia[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+
        indata$phi5[i]+indata$phi6[i]+indata$phi7[i]+indata$phi8[i]+
        indata$phi9[i]+indata$phi10[i]+indata$phi11[i]+indata$phi12[i]+
        indata$phi13[i]+indata$phi14[i]
      phib[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+ 
        indata$phi5[i]+indata$phi6[i]+indata$phi7[i]+indata$phi8[i]+
        indata$phi9[i]+indata$phi10[i]+indata$phi11[i]+indata$phi12[i]+
        indata$phi13[i]+indata$phi14[i]+indata$phi15[i]
      if(indata$phi15[i] >=0.001 & abs(1-phib[i]) > 0.001){
        theta[i] = (log(1-phia[i])-log(1-phib[i]))/(log(b)-log(a))
        k[i] = ((phib[i] - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi15[i] >=0.001 & abs(1-phib[i]) < 0.001){
        theta[i] = (log(1-phia[i])-log(1-0.999))/(log(b)-log(a))
        k[i] = ((0.999 - phia[i])/((1/a^theta[i])-(1/b^theta[i])))^(1/theta[i])
        fp[i] = 1-(k[i]/inc_lim[i])^theta[i]
      } else if(indata$phi15[i] < 0.001){
        fp[i] = phia[i]
      }
    } else if (inc_lim[i] >= 200000){
      fp[i] = indata$phi1[i]+ indata$phi2[i]+ indata$phi3[i]+ indata$phi4[i]+ 
        indata$phi5[i]+indata$phi6[i]+indata$phi7[i]+indata$phi8[i]+
        indata$phi9[i]+indata$phi10[i]+indata$phi11[i]+indata$phi12[i]+
        indata$phi13[i]+indata$phi14[i]+indata$phi15[i]#+
      #indata$phi16[i]
    } else {
      fp[i] = 0
    }
    
  }
  
  res <- cbind(geoid, cbsa, tot, fp)
  
  return(res)
  
}
